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ABSTRACT 

The theory for the formation of the first population of stars (Pop III) predicts an initial mass function (IMF) 
dominated by high-mass stars, in contrast to the present-day IMF, which tends to yield mostly stars with masses 
less than 1 M©. The leading theory for the transition in the characteristic stellar mass predicts that the cause is 
the extra cooling provided by increasing metallicity. In particular, dust can overtake H 2 as the leading coolant 
at very high densities. The aim of this work is to determine the influence of dust cooling on the fragmentation 
of very low metallicity gas. To investigate this, we make use of high-resolution hydrodynamic simulations 
with sink particles to replace contracting protostars, and analyze the collapse and further fragmentation of star- 
forming clouds. We follow the thermodynamic response of the gas by solving the full thermal energy equation, 
and also track the behavior of the dust temperature and the chemical evolution of the gas. We model four clouds 
with different metallicities (10~ 4 , 10" 5 , 10" 6 Z©, and ), and determine the properties of each cloud at the point 
at which it undergoes gravitational fragmentation. We find evidence for fragmentation in all four cases, and 
hence conclude that there is no critical metallicity below which fragmentation is impossible. Nevertheless, 
there is a clear change in the behavior of the clouds at Z = 10" 5 Z Q , caused by the fact that at this metallicity, 
fragmentation takes longer to occur than accretion, leading to a flat mass function at lower metallicities. 
Subject headings: early universe — hydrodynamics — methods: numerical — stars: formation — stars: lumi- 
nosity function, mass function 



1. INTRODUCTION 

The first burst of star formation in the Universe was thought 
to give rise to massive stars, the so-called Population III 
(Pop. Ill), with numeri cal simulations predicting masses in the 
range 20-150 M Q (e.g. [Abel et al.| [20021 |Bromm et"aL||20Q2 



|Q , Shea& Nor man, 2007 ; Yoshida et al. \ 2008). However, re 
cent results show that lower mass stars can also be fo rmed, 
albeit with characteristic masses above the solar value (Clark 



etaLj 2011abi Greif et al, 201 1] |Stacy et al. 2010 |20lT 



Smith et al.| |201 la[ |Greil et al.| |2012j). This contrasts~witn 



Tills 



present-day star formation, which ty pically yields sta rs with 
masses less than 1 M (Kroupa, 2002; Chabrier, 2003), and so 
at some point in the evolution or the Universe there must have 



) point i 

been a transition from primordial (Pop. Ill) star formation to 
the mode of star formation we see today (Pop. II/I). 

When gas collapses to form stars, gravitational energy is 
transformed into thermal energy and unless this can be dis- 
sipated in some fashion, the collapse will come to a halt. 
Thermal energy can be removed by processes such as atomic 
fine structure line emission, molecular rotational or vibra- 
tional line emission, or thermal emission from dust grains. 
In some cases, these processes are able to cool the gas sig- 
nificantly during the collapse. Th is temperature drop can pro 



|& Normanl [20071 |Yoshida~et"aLl [20081 ISmith et al.[ |201 la) . 
This happens because H2 cooling becomes inefficient for tem- 
peratures below 200K and densities above 10 4 cm~ 3 . At this 
temperature and density, the mean Jeans mass at cloud frag- 
mentation is 1 ,000 times larger than in present-day molecular 
clouds. 

Metal line cooling and dust cooling are effective at lower 
temperatures and larger densities, and so it has been proposed 
that metal enrichment of the interstellar medium by previ- 
ous generations of stars causes the transition from Pop. Ill to 
Pop. II. This suggests that there might be a critical metallicity 
Z cr i t at which the mode of star formation changes. 

The main coolants that have been studied in the l itera- 
ture are C 11 and O 1 fine structure emission (Bromm 
200H |Bromm &Loebl|2 003; Sa ntoro & Sfa ull, 2006; 



2007; 




et al.l |20071 |Jappsen et aL| 

Smi th et al. 1 120091), and dust emission (e.g. [S chneider 



2009a b ; [Smith & Sigurd sson| 



mote gravitational fra gmentation (M ac Low & Klessen 2004 
|Bonnell et al.[ |2007| ) by decreasing the Jeans mass, which 



means that instead of forming very massive clumps, with frag- 
ment masses corresponding to the initial Jeans mass in the 
cloud, it can instead form a large number of fragments with 
lower masses. 

If the gas is cooled only by molecular hydrogen emis- 
sion, numerical simulations sho w that most of the stellar mass 
would be in massive st ars (Nakam ura & Umemura[ 1 1 999 , 
|200T1 [20021 |Abel et al.l [20021 |Bromm et al.l |2002l |Q'Shea| 



etan|2002, 2006; Schn eider et aLl[20T2l[Qmukai et al.H "2005 [ 
2010). Carbon and oxygen are identified as the key species 
because in the temperature and density conditions that char- 
acterize the early phases of Pop. Ill star formation, the O i 
and C ii fin e-structure lines dominate ove r all other metal line 
transitions ( [Hollenbach & McKee]|19891 ). By equating the C 
ii or O i fine structure cooling rate to the compressional heat- 
ing rate due to free-fall collapse, one can define critical abun- 
dances [C/H ] = -3.5 and [Q/H] = -3.(Qfor efficient metal 
line cooling ( Bromm & Loeb( 2003 ). 

If one assumes that effective fine- structure cooling is a 
necessary condition for the formation of Population II stars, 
then all such stars should have a "transition discriminant" 
Arans = log(10 [C/H] + 0.3 x 10 [Q/H] ) great er than a critical 
value Arans,crit = -3.5 (Fre bel et al.[ |2QQ7] ). Although most 
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1 [X/Y] = logio(N x /N Y )* -logio(N x /N Y )o, for elements X and Y, where 
* denotes the gas in question, and where Nx and Ny are the mass fractions 
of the elements X and Y. 
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metal-poor stars lie above this value, at least one sta r has been 
obs erved to lie below it (SDSS Jl 029 15+1 72927 ; |Caffau et| 
201 lj , and there are other objects that might also have 
ues of Arans, below the c ritical value: CS30336-049 ( |Lai| 



al. 



va 



et al., 2008) and Scl07-50 (Tafelmeyer et al, 2010). 

Previous works (Jappsen et al., 2009a y have shown that 
this metallicity threshold does not represent a critical metal- 
licity. The fact that the metal-line cooling rate has a larger 
value than the compressional heating rate does not necessar- 
ily lead to fragmentation, and even in cases where it does, 
the fragments that form have masses M » 1M . This points 
towards a diffe rent process leading to low-mass stars for the 



early Universe (Kl essen et al.| 2012J ). 

A more promising way to form low mass Pop. II stars in- 
volves dust coolin g. Dust cooling models (e.g.|Om ukai et al. , 



2005 J |2QTQ1 [D opcke eTaM |20TT] [Schneide r ^^ 2006 ; 



Schneide r et al.| |2012| ) predict a much lower critical metal 



licity (Z cr it ~ 10" 4 - 10" 6 Z©), with most of the uncertainty 
coming from the nature of the dust in high-redshift galaxies. 

At densities n > 10 11 c m" 3 dust cooling becomes effi- 
cient (Omukai et al.| 2010), since inelastic g as-grain colli- 
sions are more frequent ( |Hollenbach & McKee| [T9 79 ). This 
cooling enhances fragmentation, and since it occurs at high 
densities, the distances between fragments can be very small 
dOmukail [20001 IQmukai et al.| [20051 |Schneider efaL] [2002] 
2006; Schneider & Omukai 2010| ). In this regime, interac- 
tions between fragments will be common, and analytic mod- 
els of fragmentation are unable to predict the mass distribution 
of the fragments. A full 3D numerical treatment, following 
the fragments, is needed. 

Initial attempts at modeling fragmentation in low metallic- 
ity gas were mad e by [Tsuribe & Omukai] ( |2006| |2008| ) and 
Clar k etal.| ( |20081 ). These studies described the thermal evolu- 
tion of the gas using effective equations of state derived from 
the one-zone calculations of Omukai et al. (2005), and showed 
that the cooling provided by dust does indeed lead to fragmen- 
tation. This treatment assumes, however, that the gas temper- 
ature adjusts instantaneously to a new equilibrium whenever 
the density changes and hence ignores thermal inertia effects. 
This may yield too much fragmentation. 

In Dopcke et al. ( 2011 ), we improved upon these previous 
treatments by solving the full thermal energy equation, and 
calculating the dust temperature through the energy equilib- 
rium equation. We assumed that the only significant exter- 
nal heat source is the cosmic microwave background (CMB), 
and included its effects in the calculation of the dust temper- 
ature. We found that model clouds with metallicities as low 
as 10" 4 Z« or 10" 5 Z« do indeed show evidence for dust cool- 



ing and fragmentation, supporting the predictions of Tsuribe 
& Omukai (2006, 2008) and Clark et al. (2008). 

In this work, we simulate the evolution of star-forming 
clouds for a wider range of metallicities (10~ 4 , 10" 5 , 10" 6 Z , 
and 0), and study the effect that this has on the mass function 
of the fragments that form. We also investigate how proper- 
ties such as co oling and heating rates, and number of Bonnor 



Ebert masses (Bonnor 



19561 |Ebert| |1955| > of the fragmenting 
clouds vary with metallicity and whether there is any system- 
atic change in behavior with increasing metallicity. 

2. SIMULATIONS 

2.1. Numerical method 

We model the collapse of a low-metallicity gas cloud using 
a modified version of the Gadget 2 (Springel, 2005) smoothed 



particle hydrodynamics (SPH) code. To enable us to con- 
tinue our simulation beyond the formation of the first very 
high density protostellar core, we use a sink particle approach 
(Bat e et at][T99 5 ; Jap psen et al. |2005| ), in the same way as 
in Dopck e et al.| ( |201 Tj ). Sink particles are created once the 
SPH particles are oound, collapsing, and within an accretion 
radius, h acc , which we take to be 1.0 AU. The threshold num- 
ber density for sink particle creation is 5.0 x 10 13 cm~ 3 . At the 
threshold density, the Jeans length at the minimum temper- 
ature reached by the gas is approximately one AU, while at 
higher densities the gas becomes optically thick and begins to 
heat up. Further fragmentation on scales smaller than the sink 
particle scale is therefore unlikely to occur. For further dis- 
cussion of the details of our sink particle treatment, we refer 
the reader to |Clarket aT] { |2011b) . 

We assume that the mean dust grain cross section is the 
same as for Milky Way dust and that the number density of 
dust grain s is a factor Z/ Z s maller than the Milky Way 
value (see Dopcke et al. 2011 ). To treat the chemistry and 
thermal balance of the gas, we use the same approach as in 
Clark et aL] ( |2011b]) , with the inclusion of dust cooling. The 
Clar k et al.j ( 201 lb| chemical network and cooling function 
were designed for treating primordial gas and do not include 
the chemistry of metals such as carbon or oxygen, or the ef- 
fects of cooling from these atoms, or molecules containing 
them such as CO or H2O. We justify this approximation by 
noting th at previous studies of very low-metallicity gas (e.g. 
Omu kai et al. 2005 2010| ) find that gas-phase metals have 
little in fluence on the thermal state of the gas. Omukai lk al.| 
2010) showed that H 2 and OH are efficient coolants at 



8 < n < 10 iU cm~ 3 for their one-zone model. In their hydro- 
dynamical calculations, however, the collapse is faster, and 
the effect of H 2 and OH is not perceptible. Therefore we 
do not expect oxygen-bearing molecules to have a noticeable 
effect on the thermal evolution of the gas. 

For the metallicities and dust-to-gas ratios considered in 
this study, the dominant sources of cooling are the standard 
primordial coolants (H 2 bound-bound emission and collision- 
induced emission) and energy transfer from the gas to the 
dust. Collisions between gas particles and dust grains can 
transfer energy from the gas to the dust (if the gas temper- 
ature T is greater than the dust temperature T & ), or from the 
dust to the gas (if T gY > T). Full details of th e dust cooling 
treatment can be found in Dopck e et al. ( 2011 ). 



2.2. Setup and Initial conditions 

We performed a set of four simulations, with metallicities 
Z/Z© = 10" 4 , 10" 5 , 10" 6 , and the metal-free case. Each sim- 
ulation used 40 million SPH particles. We used these simula- 
tions to model the collapse of an initially uniform gas cloud 
with an initial number density of 10 5 cm" 3 and an initial tem- 
perature of 300 K. The cloud mass was 1000 M . We in- 
cluded small amounts of turbulent and rotational energy, with 
£turb/|£gravl = 0.1 and/3 = E wt /\E gmY \ = 0.02, where £g rav is 
the gravitational potential energy, £ tur b is the turbulent kinetic 
energy, and E YOt is the rotational energy. The mass resolution 
is 2.5 x 10" 3 M o , whic h corresponds to 100 times the SPH 
particle mass (see e.g. Bate & Burkert, 1997). The redshift 
chosen was z = 15, when the cosmic microwave background 
temperature was 4 3. 6K. The dust properties were taken from 
Goldsmith (2001), and the du st grain opacities wer e calcu- 



lated in the same fashion as in |Banerjee e t al. (2006). In the 
calculations, the opacities vary linearly with Z, which means 
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Fig. l.—: Dependence of gas and dust temperatures on gas 
density for metallicities 10" 4 , 10" 5 , and 10" 6 and zero times 
the solar value, calculated just before the first sink particle was 
formed (see Table [TJ. In red, we show the gas temperature, 
and in blue the dust temperature. The dashed lines are lines 
of constant Jeans mass. 

for instance that for the Z/Z© = 10" 4 calculations, the opaci- 
ties were 10" 4 times the original values. 

3. ANALYSIS 

3.1. Thermodynamical evolution of gas and dust 

Dust cooling is a consequence of inelastic gas-grain colli- 
sions, and thus the energy transfer from gas to dust vanishes 
when they have the same temperature. We therefore expect 
the cooling to cease when the dust reaches the gas temper- 



ature. In order to evaluate the effect of dust on the thermo- 
dynamic evolution of the gas and verify this assumption, we 
plot in Figure [T] the temperature and density for the various 
metallicities tested. We compare the evolution of the dust and 
gas temperatures in the simulations, at the point of time just 
before the formation of the first sink particle (see Table [TJ. 
The dust temperature (shown in blue) varies from the CMB 
temperature in the low density region to the gas temperature 
(shown in red) at much higher densities. 

Changes in metallicity influence the density at which dust 
cooling becomes efficient. For the Z = 10" 4 Z© case, dust 
cooling begins to be efficient at n « 10 n cm" 3 , while for Z 
= 10" 5 Z©, the density where dust cooling becomes efficient 
increases to n « 10 13 cm~ 3 . For the Z = 10" 6 Z© case, dust 
cooling becomes important for n > 3 x 10 13 cm -3 , preventing 
the gas temperature from exceeding 1500 K. For comparison, 
in the metal-free case the gas reaches temperatures of approx- 
imately 2000 K. 

The efficiency of the cooling is also expressed in the tem- 
perature drop at high densities. The gas temperature decreases 
to roughly 400 K in the 10" 5 Z© simulation, and 200 K in the 
Z = 10" 4 Z case. This temperature drop significantly in- 
creases the number of Jeans masses present in the collapsing 
region, making the gas unstable to fragmentation. The dust 
and the gas temperatures couple for high densities, when the 
compressional heating starts to dominate again over the dust 
cooling. The subsequent evolution of the gas is close to adia- 
batic. 



When we compare our results to the calculations of Omukai 
et al. ( 2010J ), we find good agreement with their ID hydro- 
dynamical models, although we expect some s mall difference 
due to effects of the turbulence and rotation (see Dopcke et al. , 
201 \) and also due to the use of different dust opacity models. 



3.2. Heating and cooling rates 

The thermal evolution of the gas during the collapse takes 
different paths depending on the metallicity, as shown in the 
density-temperature diagram (Figure [TJ. In order to explain 
this behavior, we take a closer look at the cooling and heating 
processes involved. In Figure|2]we show the main cooling and 
heating rates divided into four panels for the different metal- 
licities. These rates were calculated by averaging values of 
individual SPH particles in one density bin, where the total 
density range was divided in 500 bins in log space. 

At densities below n « 10 10 cm~ 3 , dust cooling is unimpor- 
tant in all of the runs. At these densities, the dominant coolant 
is H2 line emission, while the heating is dominated by com- 
pressional (pdV) heating at n < 10 8 cm~ 3 , and by three body 
H2 formation heating at higher densities. 

At higher densities, dust cooling starts to play a more im- 
portant role. In the Z = 10" 4 Z o simulation, dust cooling ex- 
ceeds pdV heating at n « 10 10 cm" 3 , although it does not 
exceed the H2 formation heating rate until n « 10 n cm" 3 . 
Once this occurs, and dust cooling dominates, the gas tem- 
perature drops sharply. In the Z = 1O" 5 Z simulation, on 
the other hand, dust becomes the dominant coolant only at 
n « 10 13 cm~ 3 , and so the temperature decrease happens later 
and is smaller. Finally, in the Z = 10" 6 Z o case, dust cooling 
becomes competitive with pdV heating only at the very end 
of the simulation, and so the effect on temperature evolution 
is less evident. 

The other thermal processes play a minor role during the 
collapse. For example, H2 dissociation cooling only becomes 
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important in the runs with Z = 10 6 Z© and 0, and only for 



10 13 cm" 3 



n > l(J 1J cm~ J . At very high densities (n > 10 cm" 3 ), H2 
collision-induced emission (CIE) cooling also begins to be 
important. For more details on H2 heating and cooling pro- 
cesses in this very high density regime, we refer to |Clark et| 
|aEl (2011b) . 

3.3. Fragmentation 

The transport of angular momentum to smaller scales dur- 
ing the collapse leads to the formation of a dense disk-like 
structure, supported by rotation. This disk then fragments into 
multiple objects. 

Figure [3] shows the density structure of the gas immediately 
before the formation of the first protostar. The top-left panel 
shows a density slice on a scale comparable to the size of 
the initial gas distribution. The structure is very filamentary 
and there are two main over-dense clumps in the center. If we 
zoom in on one of the clumps, we see that its internal structure 
is also filamentary. Observe that at large scales the gas cloud 
properties are the same for all metallicities. Differences in the 
thermodynamic evolution appear only at n > 10 11 cm -3 (see 
Figure Q. As a consequence, we observe variations in the 
cloud structure only in the high-density regions. 

Once the conditions for sink particle creation are met (see 



Section 
(Figure 



2.1|), they start to form in the highest density regions 
Then, a disk is built up in the se regions, where 



1980). During further 



fragmentation also occurs (Tohline 
collapse, this dense region creates spiral structures. For Z 
= 10" 5 Z Q and 10" 4 Z Q , density waves build up spiral struc- 
tures, which become locally gravitationally unstable and go 
into colla pse. The formation of bina ry systems by triple en- 
counters (Binn ey & Tremaine} |2008| ) transfers kinetic energy 
to some sink particles, causing them to be ejected from the 
high density region. For Z = 10" 5 Z Q , when the star formation 
efficiency (SFE) is 0.5%, fragmentation has already occurred 
in a secondary dense center, at a distance of ~ 20 AU from the 
first dense region. 

For Z = 10" 6 Z© and 0, the formation of spiral structures is 
not observed. In these two runs, star formation occurs mainly 
in the central clump. 

One way to study the effect of dust cooling on the frag- 
mentation behavior and the final stellar IMF is to look at the 
changes in the number of Bonnor-Ebert (M BE ) masses con- 
tained in this centr al dense region. Using the definition from 
|Brommetal.l ( |2009] ), 

M - = 500M °(200k) (low) ' (1) 

for an atomic gas with temperature T and number density n, 
we have computed the number of Bonnor-Ebert masses con- 
tained within a series of concentric radial spheres centered on 
the densest point in each of our four simulations. The results 
are shown in Figure [4| 

At the beginning of the simulation, the cloud had ~ 3 M BE . 
During the collapse, the gas cools and reaches ~ 6 M BE in 
all cases. Cooling and heating are different depending on the 
metallicity, and this difference is seen for distances smaller 
than ~ 400 AU. The Z = 10" 4 Z© case, for instance, has 
twice the number of M B e for distances smaller than ~ 10 AU, 
when compared to the other cases. This will have direct con- 
sequences for the fragment mass function as we will see in the 
next section. 
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TABLE 1 : Sink particle properties for the different metallici- 
ties at the point where 4.7 M© have been accreted by the sink 
particles. "ST" (start time) is the time when sink particles 
start to form. "FT" (formation time), is the time taken to ac- 
crete 4.7 M Q in the sinks. "SFR" is the mean star formation 
rate. Mean and median refer to the final mean and median sink 
mass. Finally, "N" is the number of sink particles formed. 

3.4. Properties of the fragments 

The simulations were stopped at a point when 4.7 M© of 
gas has been accreted into the sink particles, because the high 
computational cost made it impractical to continue. Figure 
[6] shows the mass distribution of sink particles at that time. 
We typically find sink masses below 1 M , with somewhat 
smaller values in the 10" 4 Z© case compared to the other 
cases. No sharp transition in fragmentation behavior was 
found, but rather a smooth and complex interaction between 
kinematic and thermodynamic properties of the cloud. 

Table [T] lists the main sink particle properties. It shows that 
the time taken to form the first sink particle is slightly shorter 
for higher metallicities. This shorter time is a consequence of 
the more efficient cooling by dust, which decreases the ther- 
mal energy that was delaying the gravitational collapse. In 
Table Q] we also observe that the star formation rate is lower 
for Z = 10" 4 Z . This is because star formation started at 
an earlier stage of the collapse, when the mean density of the 
cloud was lower and there was less dense gas available to form 
stars. 

To better understand whether the resulting stellar cluster 
was affected by varying the metallicity, we plot the final sink 
mass distribution in Figure|6] It shows that for the simulations 
with Z < 10" 5 Z©, the resulting sink particle mass function is 
relatively flat. There are roughly equal numbers of low-mass 
and high-mass stars, implying that most of the mass is to be 
found in the high-mass objects. This mass function is consis- 
tent with those found in other recent studies of fragmenta tion 
in metal-free gas ( |Greif et al.[ |201 1[ |Smith et al.| |201 la| ). If 
the sink particle mass function provides a reliable guide to the 
form of the final stellar IMF, it suggests that at these metallic- 
ities, the IMF will be dominated by high-mass stars. 

All of the histograms in Figure |6] have the lowest sink par- 
ticle mass well above the resolution limit of 0.0025M©. Note 
that in all cases, we are still looking at the very early stages 
of star cluster evolution. As a consequence, the sink particle 
masses in Figure [6] are not the same as the final protostellar 
masses - there are many mechanisms that will affect the mass 
function, such as continuing accretion, mergers between the 
newly formed protostars, feedback from winds, jets and lumi- 
nosity accretion, etc (see Section]?]). 

3.5. Timescales 

One way to explain the final mass distribution of the frag- 
ments is to look at the timescales for mass accretion and frag- 
mentation. The degree of gravitational instability inside a 
volume can be represented by the number of Bonnor-Ebert 
masses contained in this volume. We can therefore estimate 
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Fig. 2.—: Cooling and heating rates versus number density for Z = 10" 4 , 10" 5 , 10" 6 Z Q , and zero. The values are calculated just 
before the first sink formed. The lines labeled as "C" indicate cooling, and "H" is heating. "Dust C", "H2 Line C", "H2 CIE" , 
and "H2 Diss." indicate dust grain cooling, H2 line emission, collision-induced emission, and dissociation cooling, respectively. 
"H2 Form. H" and "pdV H" are the H2 formation heating rate, and compressive (pdV) heating rate. 



the fragmentation timescale by computing the time taken for 
the central dense region to accrete one Bonnor-Ebert mass. In 
other words, we have tf mg = M B e/M, where M is the accre- 
tion rate. This value is shown as a function of the enclosed gas 
mass in Figure [7] where the values are calculated for particles 
in spherical shells, and the center is taken to be the densest 
SPH particle. M is obtained by summing up the mass of the 
particles (m p ) inside a shell times their radial velocity (v r ), 

M=2_j m P V r- 
shell 

For comparison, we also plot the accretion timescale, here 
defined as the time taken by the gas to accrete the mass en- 
closed by that radius, t acc = M enc /M. When the fraction 
feag/^acc > 1> one expects that the gas enclosed by this shell 
is going to be accreted faster than it can fragment, favoring 
high mass objects. Conversely, for tf V3Lg /t acc < 1, the gas will 



fragment faster than it can be accreted by the existing frag- 
ments, and the final mass distribution is expected to have more 
low mass objects. Note that as defined here, the timescale on 
which new fragments form tf mg /t acc is the inverse of the quan- 
tity M gas /M B E plotted in Figure [?| 

In Figure [7J the simulation with Z = 10" 4 Z© has the lowest 
values for tf mg /t acc , over a wide range of M enc . This indicates 
that more low-mass fragments are expected to form in this 
case, leading to a steeper fragment mass function. 

Now we can compare the predicted values before sink 
formation started with the final accretion and fragmentation 
timescales. These values are designed to represent the charac- 
teristic timescales on which the mass histogram changes: the 
fragmentation timescale (rf rag ) is the time on which the num- 
ber of fragments change by a significant amount, while the ac- 
cretion timescale (T acc ) represents the time on which the exist- 
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Fig. 3.—: Number density maps for a slice through the high 
density region for Z = 10" 4 Z© (top), 10" 5 Z , 10" 6 Z , and 
(bottom). The image shows a sequence of zooms in on the 
density structure in the gas immediately before the formation 
of the first protostar. 
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Fig. 4.—: Enclosed gas mass divided by Bonnor-Ebert mass 
versus radius for different metallicities. The values were cal- 
culated at the time just before the first sink was formed and the 
center is taken to be the position of the densest SPH particle. 



ing fragments grow in mass. We therefore define Tf rag = n/h, 
and r acc = M/M, where n is the number of sink particles, and 
M is the total mass incorporated into sink particles. 

Both timescales should increase over time, since the num- 
ber of sink particles and the total mass also increase. This 
reflects the fact that it takes longer to change the shape of 
the mass histogram in frequency when the number of ele- 
ments is higher, and in mass when the total mass is high. The 
key point is that these timescales are related to the shape of 
the histogram. If the timescale for accretion is shorter than 
the timescale for fragmentation, the histogram will tend to be 
dominated by high-mass objects. Conversely, if the fragmen- 
tation timescale is shorter than the accretion timescale, the 
low-mass part will be populated before the objects can grow 
and occupy higher mass bins. Therefore, a comparison be- 
tween r aC c and Tf r ag helps us to understand the shape of the 
mass histogram. 

We calculate the average r acc and Tf rag at each point in time 
when a new sink particle is created, since Tf rag is not defined 
for h = 0. The value for the mean Tf rag and r acc is calculated 
by averaging their individual times in equations [2] and |3j 



Tfrag,!-,.. - 
7"acc,i — • — 



Ant/Ati 

Mt 



n t - fli-i 
Mi{ti-ti.i) 



AMi/Ati Mi - M, 



(2) 
(3) 



i-\ 



where / is the point in time when a new sink is created, and it 
varies from 2 to the total number of sink particles, n^ M t , tt 
are the number of sink particles, the total mass in sinks, and 
the time at the point i. 

Figure [8] shows the average timescales for fragmentation 
and accretion for different metallicities. These results explain 
the difference in the sink particle mass distribution in Figure 
[6] For Z < 10" 5 Z©, the fragmentation time is always larger 
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Fig. 5.—: Number density map showing a slice through the densest clump, and the star formation efficiency (SFE) for Z = 10" 4 Z Q 
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3.— : Sink particle mass function at the point when 4.7 
M Q of gas had been accreted by the sink particles. The mass 
resolution of the simulations is indicated by the vertical line. 



than the accretion time, indicating that the sink particles will 
accrete faster than they can be generated, resulting in a flatter 
mass distribution. On the other hand, when the fragmenta- 
tion time is longer than the accretion time (for Z = 10" 4 Z ), 
the gas fragments, rather than moving to the center and being 
accreted. As a consequence, the low-mass end of the proto- 
stellar mass function grows faster than the high-mass end, and 
the slope of the mass function steepens. This behavior agrees 
well with the predictions from before fragmentation started, 
shown in Figure [7] 

Note that the values in Figure [7] were calculated before the 
formation of the first sink particle, while the values in Figure 
[8] were calculated using the sink particle properties. A com- 
parison between them is useful to evaluate whether the gas 
cloud properties from before star formation started could be 
used to predict its star formation behavior. The trend in both 
figures is that the fragmentation timescale is normally shorter 
than the accretion timescales for Z = 10" 4 Z . From these 
results we conclude that the mass distribution in Figure [6] can 
be explained by the timescales in Figure [8] in particular the 
fact that the Z < 10" 5 Z© simulations have more high-mass 
objects. The last finding is that the transition from Pop. Ill 
to Pop. II star formation mode is not abrupt, in the sense that 
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Fig. 7.—: Timescales for fragmentation (top panel) and ac- 
cretion (middle panel), and also their ratio (bottom panel) ac- 
cording to the radius for the metallicities tested. The values 
were calculated just before the first sink particle was formed. 
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Fig. 8.—: Timescales for fragmentation and accretion for dif- 
ferent metallicities calculated for the sink particles following 
Equations [2] and [3] 

there is no metallicity bellow which the gas cannot fragment. 
The transition is rather in the stellar initial mass function, and 
gas clouds with Z < 10" 5 Z© form a more flat IMF, while gas 
clouds with Z > 10~ 4 Z produce a cluster with more low- 
mass objects (see also Clark et aLl |2008). 



3.6. Radial mass distribution 
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Fig. 9.—: Dependence of the enclosed gas and sink mass on 
the distance from the center of mass, for the four simulations. 
The values were calculated at a point when 4.7 M© of gas had 
been accreted. 

Another property of the star-forming cloud that we ob- 
served to vary in our calculations is the spatial mass distri- 
bution. The dependence of the enclosed gas and sink mass 
on the distance from the center of mass is shown in Figure 
[9] The Z = case has almost all of the sink particle mass 
concentrated within r < 8AU. The gas density for this case is 
also higher in this region, when compared to the other metal- 
licities, showing that the gas and sink particle mass densities 
follow each other. The mass in sink particles exceeds the gas 
mass for small radii, being the most important component in 
the gravitational potential. For r > 15 AU, the gas becomes 
the most massive component, for all Z. |Girichidis etal. ( 2012 ) 
also reported this behavior, but in their case the sink particles 
already started to dominate the potential below r « 10 3 AU. 

This higher concentration of gas and sinks at the center oc- 
curs because for the Z = case, the gas had higher tempera- 
tures in the central region. For high temperatures, the criterion 
for gravitational instability requires higher densities, which 
are achieved only very close to the center. As a consequence, 
the sink particle formation criteria are met just for short dis- 
tances from the center. 

Consequently, the dominance of sink particles mass in the 
gravitational potential over the gas mass, for radii smaller than 
150 AU, shows the importance of treating gas and stars to- 
gether in this sort of problem. It also suggests that N-body 
effects, such as ejections and close encounters, should play an 
important role in the formation of these dense star clusters, 
eve n in the very earliest stage s of their evolution (see |Smithet 
|aEl [20TTH |Greif et al.| |20l2) . 

4. CAVEATS 

Our aim with these calculations is to study the importance 
of dust cooling for fragmentation in high-redshift halos. To 
better understand star formation in this environment, addi- 
tional physical processes should be considered as well. 

Particularly, the low number of sink particles (« 20) and the 



small SFE (0.5%) do not permit to constrain the stellar IMF 
adequately. By running the calculations until the SFE goes 
towards higher values, uncertainties involving the fragments 
that formed during the simulations can be diminished. How- 
ever, this does not appear to be computationally feasible with 
our current approach. 

The mass accreted by the sink particles varies with the dif- 
ferent metallicities, and affects the final sink particle mass 
function. It also influences the expected accretion luminosity. 
We did not take this process into account during our calcu- 
lations, but we can estimate its importance relative to other 
thermal processes. 
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Fig. 10.—: Time evolution of the mass, mass accretion rate, 
and accretion luminosity for the four metallicities, for all sink 
particles combined. 



In Figure 
newborn ste 



10] we present the accretion properties for the 
Tar systems. The top panel shows how the to- 
tal mass in sinks evolve with time, for different metallicities. 
The accretion rate varies from 0.02 to 0.17 M© yr -1 , and it is 
on average lower for the Z = 10" 4 Z n case. 



In the bottom panel of Figure 10 we show the accretion 
luminosity calculated by adding up all sink particle contribu- 
tions, with the standard equation, 



GM*M 



(4) 
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where M is the mass accretion rate by a protostar wi th mass 
M*, and s tellar radius R*. We calculate R* following Stahler 



etal.| { |1986) using 



R* = 66.8 






0.27 



M 



lO" 2 M yr- 



0.41 



R*. 



(5) 



The accretion luminosity varies from few times 10 3 L Q to 
around 50 x 10 3 L Q , depending on metallicity. For Z < 10" 5 
Z , the accretion luminosity is always over 10 4 L Q . The Z 
= 10" 4 Z case has the lowest estimated accretion luminos- 
ity, around four times lower than the other cases. The val- 
ues found for the accretion l uminosity are similar to the ones 
found by Smit h et al. ( 20 1 1 a| ) for Z = 0, where they argue that 
accretion luminosity could delay the fragmentation, but not 
prevent it. 

We can now compare this expected luminosity, and its con- 
sequent heating, to the heating processes in Figure [2] We 
make the assumption that gas and dust are absorbing the ra- 
diation in the optically thin regime. This overestimates the 
effects, and we obtain an upper limit for the accretion lumi- 
nosity heating, 



T^ / L* acc \ i 



(6) 



where p g is the gas density, Kp is the Planck mean opacity, 
and r is the distance from the source. We also assume that 
heating occurs atp « 10" 9 g cm -3 . With this assumption we 
can calculate the mean gas temperature, and use the Planck 
mean opacity for the gas from Mayer & Duschl (2005), for 
their fiducial Pop. Ill chemical composition. The Plan ck mean 
opacity for th e dust is calculated in the same way as in Dopcke 
et al. ( 2011| ). Finally, the combined Planck mean opacity is 
the sum of gas and dust contributions, k p = K gas + Kd U st- 

By considering the maximum accretion luminosity for 
each case, we get that Y acc = 4.9, 0.9, 1.7, and 0.7 x 
10 3 (20AU/r) 2 erg s^g -1 , forZ = 10" 4 , 10" 5 , 10" 6 , andOZ . 

As these values are comparable to the other thermal pro- 
cesses at high densities (see Figure [2]), it would seem that 
accretion luminosity heating from the young protostars may 
have some effect on the way in the which the gas behaves. 
However without doing the radiative transfer explicitly, it is 
difficult to estimate how big this effect will be. 

Although the amount of heating seems high, the dust cool- 
ing is a strong function of temperature in this regime, and so 
it could be that dust temperatures remain quite similar. One 
must also remember that the above estimates do not take into 
account the extinction and reprocessing of the radiation field 
that will occur in the optically thick region that surrounds the 
protostar. However even a factor of 2 change in the dust tem- 
perature will remove the dip in the p - T phase diagrams that 
we show in Figure [T] and thus remove the ability of the dust 
to set a new length-scale for fragmentation. 

Sufficiently far enough away from the strongest sources, the 
effect will obviously drop to the point at which the physics in 
our current calculations are applicable. However if we look to 
Figure[5] we see that most of the fragmentation that we report 
is confined to a few tens of AU around the central protostar, 
and as such, the effects of the accretion luminosity are likely 
to change the picture that we present in this paper to some 
extent. We hope to explore this effect in a future study. 

Another aspect of our model that could be improved upon 
are the dust opacities. The thermal evolution can be calculated 
more accurately if we use dust opacities that better represent 



the values expected for very low metallicity environments. 
The dust opacity in our simulations correspond to values cal- 
culated for the Milky Way and then scaled with metallicity. 
This means that the opacity values for the Z = 1O" 4 Z case are 
10" 4 times the dust opacity in the Milky Way. This approxi- 
mation is probabl y not f ully correct, and the use of a mo re ac- 
curate model (e.g. Todini & Ferrara, 2001 ; Bianchi & Schnei- 
[der[ [20071 jSchneider et al.||2012[|Nozawa et 31^2(305^ 2006) 
can change the value of the cooling in the region where frag- 
mentation occurs. This change affects the local Jeans mass, 
and consequently the star formation behavior. 

Furthermore, the available models give the dust composi- 
tion for different scenarios in the early Universe, e.g. dif ferent 
supernovae progenitor masses (Schneide r et al.| |2012[ ), and 
the use of such models would add another variable to the prob- 
lem - the stellar population for the supernovae progenitors. 
One reasonable approach is to test different scenarios and see 
how they would affect the properties of the cluster of stars 
that forms. In this sense, the dust composition is a problem in 
itself that should be addressed. Since cooling affects the frag- 
mentation behavior and mass accretion, a more realistic dust 
model improves the accuracy with which we can model star 
formation at low metallicities. We intend to address this issue 
in a future paper. 

Another effect that could change our results is the possibil- 
ity of inelastic encounters between the protostars. Star for- 
mation in our simulation occurs at very high densities, where 
inelastic encounters between the new born protostars could 
occ ur. In sim ilar conditions to the ones tested here, |Smith et| 
|aL| ( |2011a"|b] ) show that the estimated stellar radius could be as 
large as ~ 1 AU, a value comparable to the distances between 
the sink particles shown in Figure [5] By not accounting for 
merging of such objects, we could be overestimating the final 
number of fragments, although we expect new protostars to 
conti nue to form. For a detailed investigation of this process 
see Greifetal.| ( |20T2l ). 

Finally, the inclusion of magnetic fields in the calculations 
could alter the fragmentation picture as it is presented in this 
study. They can be am plified during gravitational collapse 
(Schleicher et al., 2010), generating values strong enough to 
delay the collap se jSchleicher eTat) 2009; S ur et al.| [20l0| 
Federrath etaE) |2011a|bj|Turk et al ||2012|). Analytic ampli- 



fication values are calculated by Schober et al. ( 2012b . From 
modeling present-day star formation, we know that the pres- 
ence of magnetic fields can decrease the level of fragmenta- 
tion, but cannot prevent it, for the expected saturation level s 
of a few percent ( Peters et aL) 201 l[ |Hennebel le et al. 2011| ). 



5. CONCLUSIONS 

In this paper we have addressed the question of whether 
dust cooling can lead to the fragmentation of low-metallicity 
star- forming clouds. For this purpose we performed numeri- 
cal simulations that follow the thermodynamical and chemical 
evolution of collapsing clouds. The chemical model included 
a primordial chemical network together with a description of 
dust effects, where the dust temperature was calculated by 
solving self-consistently the thermal energy equilibrium equa- 
tion. 

As a result, we found that dust can cool the gas, for num- 
ber densities higher than 10 11 , 10 12 , and 3 x 10 13 cm~ 3 for Z 
= 10" 4 , 10" 5 , and 10" 6 Z , respectively. Higher metallicity 
implies larger dust-to-gas fraction, and consequently stronger 
cooling. This is reflected in a lower temperature of the dense 
gas for the higher metallicity simulations, and this colder gas 
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permitted a faster collapse. Therefore, the fragmentation be- 
havior of the gas depends on the metallicity, and higher metal- 
licities lead to a faster collapse. 

For example, the characteristic fragment mass was lower 
for Z = 10" 4 Z , since a lower temperature reduces the 
Bonnor-Ebert masses at the point where the gas undergoes 
fragmentation. This also implies a lower ratio of fragmenta- 
tion and accretion time, tf mg /t acc , which will lead to a mass 
function dominated by low-mass objects. For Z < 10" 5 Z , 
fragmentation and accretion timescales are comparable, and 
the resulting mass spectrum is rather flat, with roughly equal 
numbers of stars in each mass bin. 

In addition to that, dust cooling appears to be insufficient to 
change the stellar mass distribution for the Z = 10" 5 and 10~ 6 
Z cases, when compared with the metal-free case. This can 
be seen in the sink particle mass function (Figure [6]), which 
shows that the Z < 10" 5 Z cases do not appear to be funda- 
mentally different. 

Finally, we conclude that the dust is not an efficient coolant 
at metallicities below or equal to Z cr i t = 1O" 5 Z , in the sense 
that it cannot change the fragmentation behavior for these 



metallicities. Our results support the idea that low mass frag- 
ments can form in the absence of metals, and clouds with Z < 
Zcrit will form a cluster with a flat IMF. 
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